Anticancer derivative of the natural alkaloid, theobromine, inhibiting EGFR protein: Computer-aided drug discovery approach

A new semisynthetic derivative of the natural alkaloid, theobromine, has been designed as a lead antiangiogenic compound targeting the EGFR protein. The designed compound is an (m-tolyl)acetamide theobromine derivative, (T-1-MTA). Molecular Docking studies have shown a great potential for T-1-MTA to bind to EGFR. MD studies (100 ns) verified the proposed binding. By MM-GBSA analysis, the exact binding with optimal energy of T-1-MTA was also identified. Then, DFT calculations were performed to identify the stability, reactivity, electrostatic potential, and total electron density of T-1-MTA. Furthermore, ADMET analysis indicated the T-1-MTA’s general likeness and safety. Accordingly, T-1-MTA has been synthesized to be examined in vitro. Intriguingly, T-1-MTA inhibited the EGFR protein with an IC50 value of 22.89 nM and demonstrated cytotoxic activities against the two cancer cell lines, A549, and HCT-116, with IC50 values of 22.49, and 24.97 μM, respectively. Interestingly, T-1-MTA’s IC50 against the normal cell lines, WI-38, was very high (55.14 μM) indicating high selectivity degrees of 2.4 and 2.2, respectively. Furthermore, the flow cytometry analysis of A549 treated with T-1-MTA showed significantly increased ratios of early apoptosis (from 0.07% to 21.24%) as well as late apoptosis (from 0.73% to 37.97%).


Introduction
The world health organization (WHO) anticipated that during the next few years, cancer will dominate all other causes of death [1]. Developing treatments that suppress the growth of cancer by interacting with specific molecular targets and damaging the cancer cells is a major  [3,4-d]pyrimidine that showed excellent efficacy for inhibiting EGFR-TK at nono-molar doses [34,35]. Our team previously synthesized compound V (a thieno [2,3-d]pyrimidine derivative) that was promising anti-proliferative and EGFR inhibitor [36] (Fig 1). These compounds possess some pharmacophoric features of EGFR-TKIs. These features are a planar heterocyclic system, an NH spacer, a terminal hydrophobic head and a hydrophobic tail. The key roles of the above-mentioned structural moieties are to occupy the adenine binding pocket [37], interact with amino acid residues in the linker region [38], to be inserted in the hydrophobic region I [39], and to occupy the hydrophobic region II [40,41], respectively (Fig 1).
In this work and as an extension of our previous efforts in the discovery of new anti-EGFR agents [36,[42][43][44], compound V was used as a lead compound to reach a more promising anticancer agent targeting EGFR. Several chemical modifications were carried out at four positions. The first position is the planar heterocyclic system. We applied the ring variation strategy as the thieno [2,3-d] pyrimidine moiety was replaced by a xanthine derivative (3-methyl-3,7-dihydro-1H-purine-2,6-dione). The six hydrogen bond (HB) acceptors may facilitate the HB interaction in the adenine binding pocket. Chain extension strategy was applied in the liker region through the replacement of the NH-linker with acetamide moiety. The terminal hydrophobic head (3-iodobenzoic acid) of the lead compound was replaced by toluene moiety) via ring variation strategy. A simplification strategy was applied for the hydrophobic tail (cyclohexene) of the lead compound. It was replaced by methyl group at 7-posision of xanthine moiety (Fig 2).
The 1 H NMR spectrum of T-1-MTA showed singlet signal at δ = 8.07 for CH imidazole and multiplet signals ranging from δ 7.41 to 6.87 for aromatic protons besides remarkable singlet signals for the CH 3 (of m-tolyl group) and CH 2 groups at δ = 2.27 and 4.67, respectively. The IR spectrum of the same product revealed absorption bands at 1711, 1662 cm -1 corresponding to carbonyl groups and absorption bands at 3255 cm -1 corresponding to NH. Regarding the 13 C NMR spectrum, four shielded signals appeared at 43.84, 33.66, 29.90, and 21.62 ppm corresponding to CH 2 and the three CH 3 groups, respectively.

Molecular docking
The examined proteins' X-ray structures (EGFR WT ; PDB: 4HJO and EGFR T790M ; PDB: 3W2O) were acquired from the Protein Data Bank (PDB, http://www.pdb.org). First, the docking protocol was verified for both wild and mutant EGFR and the RMSD results were 1.20 and 1.15 Å, respectively Erlotinib, as a native inhibitor for EGFR WT , revealed an affinity value of -20.50 kcal/mol. The binding pattern of erlotinib revealed a key HB with Met769 (2.11 A˚) in addition to four hydrophobic interactions (HI) in the adenine pocket and three HIs with Ala719 and Val702, and Lys721 in the hydrophobic pocket (Fig 4). TAK-285, as a native inhibitor for EGFR T790M , presented a binding energy of -7.20 kcal/mol. The binding pattern of TAK-285 revealed a key HB with Met793 (2.44 A˚) through the pyrimidine moiety in the adenine pocket. The later moieties (3-(trifluoromethyl)phenoxy and N-ethyl-3-hydroxy-3-methylbutanamide moieties)

PLOS ONE
were fixed in the hydrophobic pocket via a network of HIs with Lys745, Ile759, Met790, Val726, and Ala743, and Leu844 (Fig 5).
Regarding the EGFR WT , a comparable affinity value to erlotinib was obtained by T-1-MTA (-20.45 kcal/mol). Additionally, it interacts with the EGFR WT active site similar to erlotinib and adopts the same orientation. Besides, the 3,7-dimethyl-2,6-dioxo-2,3,6,7-tetrahydro-1Hpurine arm formed a crucial HB with Met769 besides two HIs with Lue694 inside the adenine pocket. On the other side, five HIs with Leu764, Ala719, Val702, and Lys721 were achieved via the m-tolyl moiety in the conserved hydrophobic pocket. The methyl group at 7-posision of xanthine moiety failed to form HIs in the hydrophobic pocket II (Fig 6).

MD simulations
The MD analyses obtained on a 100 ns production run showing an overall system stability. The RMSD plot (Fig 8A) showed a stable trend for the EGFR only and the EGFR_T-1-MTA complex that were represented as blue and green curves showings averages of 2.16 Å and 2.97 Å, respectively. Moreover, the RMSD of the T-1-MTA (red) showed three states during the whole trajectory. The first 10 ns show an average of 2.16 Å before spiking to an average of 9.43 Å for the next 30 ns. Moreover, the last 60 ns show a large stable average value of 17.72 Å. The reason for this increase in the RMSD values of the compound T-1-MTA is due to the translational movement of the compound T-1-MTA relative to the protein as shown in Fig 8G which compares between the positions of the ligand at 1.5 ns (green sticks), 29.5 ns (cyan sticks), 83.9 ns (magenta sticks), and 94 ns (yellow sticks). The RoG (Fig 8B), SASA and (Fig 8C) HB show a stable protein fluctuation with an average of 19.51 Å, and 15285 Å 2 , respectively. The change in HBs between the T-1-MTA and EGFR (Fig 8D) shows that there is, approximately, at least one HB formed during the first 40 ns and it increases to at least two bonds during the rest of the simulation. The amino acids' fluctuation was depicted in the RMSF plot (Fig 8E) showing low values of fluctuation (less than 2 Å) excepting the free C-terminal and the loop region E842:Y845 reaching 7 Å, and 3.5 Å, respectively. During the simulation time, the distance between the center of mass of compound T-1-MTA and the center of mass of EGFR protein shows a similar trend to the RMSD values of the ligand (three states) (Fig 8F). It started with an average of 16.72 Å for the first 15 ns before slightly decreasing to an average of 14.02 Å for the next 25 ns (from 15 ns to 40 ns). Finally, the last 60 ns showed an average value of 11.87 Å showing a stable interaction (Fig 8G).

MM-GBSA studies
The binding free energy of the EGFR_T-1-MTA complex was further analyzed deeply by the MM-GBSA analysis. As Fig 9 shows, the EGFR_T-1-MTA complex had a total binding energy

Protein-Ligand Interaction Profiler (PLIP) studies
After that, to obtain a representative frame for each cluster of the EGFR_T-1-MTA complex, the obtained trajectory was clustered. The elbow method was used to automatically choose the number of clusters, as described in the methodologies section, and this resulted in four clusters. The PLIP website was used to determine the number and types of interactions between T-1-MTA and EGFR for each cluster representative ( Table 1). As can be seen, HIs have a similar overall number of interactions in all the clusters compared to the HBs (7 HIs vs. 6 HBs). Additionally, a.pse file was generated to understand the 3D conformations of T-1-MTA as well as its interaction against the EGFR (Fig 11).

DFT studies
In an attempt to clarify the inhibitory activity of T-1-MTA, theoretical DFT studies have been explored. The conceptual DFT has been used for understanding the electronic structure of the prepared molecule to determine its structural features which has far-reaching consequences on the molecules' reactivity. Hence, the DFT-based reactivity descriptors (global), frontier molecular orbital analysis (FMO), and surface potential maps have been investigated to explore the reactivity of the prepared compound.

Geometry optimization
The reactivity of T-1-MTA is mainly determined by its chemical structure, so the structure is fully optimized and computed using DFT. The single bond length N2-C14 is 1.4765 Å,

Frontier molecular orbital analysis (FMO) analysis
Border molecular orbitals in a molecule play a vital role in the electric properties as the system with a smaller value of energy gap between the border orbitals (E gap = E LUMO -E HOMO ) should be more reactive than one having a greater E gap . Fortunately, T-1-MTA reported a smaller E gap value, so the electronic movement between the border orbitals; LUMO and HOMO, could occur easily [47]. The nodal properties of HOMO-LUMO orbitals of the studied heterocyclic molecule in

PLOS ONE
within the molecule is easy [48]. The quantum chemical parameters such as ionization potential (IP) and electron affinity (EA) were calculated and listed in Table 2.
Global reactive indices and total density of state (TDOS). Based on the density functional theory (DFT) concept, global reactivity parameters are essential tools for comprehending the behavior of any chemical molecular structure. Such global reactivity indices depend on the value of E gap . In Table 2, the static global properties of T-1-MTA, namely the electrophilicity (ω), maximal charge acceptance (N max ), energy change (ΔE), chemical potential (μ), global chemical softness (σ), global electronegativity (χ), global chemical hardness (η), and electron   Table 2 indicated that T-1-MTA is treated as soft within the nucleophilicity and electrophilicity scales [49]. The density of states and the distribution function probability determined by the occupied states per unit volume are important to provide an accurate description best than frontier molecular orbitals. The TDOS spectrum of T-1-MTA in Fig 14 depicted that the highest electronic intensity is located in the occupied orbitals under the HOMO orbital. Also, the TDOS spectrum confirmed the narrow HOMO-LUMO gap.
2.6.4. Molecular surface potential maps. Molecular electrostatic surface potential discovers the relationship between the electronic distribution over the molecule surface and its binding ability. The molecular electrostatic potential explains and predicts the noncovalent

ADMET profiling study
The approval of any new compound as a marketed drug is based on a pharmacokinetic evaluation in addition to its biological activity. So, analyzing the ADME properties of a compound at the early stages should keep the discovery process from being delayed [50]. Although ADMET studies in vitro can investigate the properties of the absorbent, distribution, metabolism, excretion, and toxicity of drugs, in silico studies are advantageous because of their ability of saving cost, time, effort in addition to the regulations restricting the use of animals [51]. Computing

PLOS ONE
ADMET parameters using Discovery is used to determine the ADMET parameters for T-1-MTA against erlotinib. Interestingly, the obtained results of T-1-MTA comparing erlotinib (Fig 16 and Table 3) showed a high likeness degree as it was anticipated to have a low potential to pass the BBB. Additionally, hepatotoxicity (HT) and the inhibition of cytochrome P-450 (CYP2D6-I) were expected to be absent. Also, T-1-MTA levels of aqueos solubility (AS) and intestinal absorption (IA) were computed as good.

In silico toxicity studies
For a drug to be developed successfully, toxicity assessment at the early stages must be done in order to control the possibility of failure in the clinical stage [52]. The in silico approach to toxicity assessment is promising being accurate and avoiding ethical and resource constraints in the in vitro and in vivo phases of toxicity development [53]. In silico prediction of toxicity basically uses the structure-activity relationship (SAR)-predicting toxicity. In detail, the computer compares the chemical properties of the examined molecules against the structural properties of tens of thousands of compounds of reported safety or toxicity [54]. Employing the Discovery studio software, eight toxicity models were used to estimate T-1-MTA's toxicity in  Table 4) 2.9. Biological evaluation 2.9.1. In vitro EGFR inhibition. For the purpose of examining the design and the computational outcomes that clearly demonstrated T-1-MTA's significant affinity for EGFR, T-1-MTA's inhibitory ability was assessed in vitro against the EGFR protein (Fig 17). The obtained inhibition value (22.89 nM) was near to erlotinib's value, and the resulting in vitro results confirmed T-1-MTA 's suppressive potential.

Cytotoxicity and safety
In vitro cytotoxicity assessment was performed for T-1-MTA using compared to erlotinib as demonstrated in Table 5. The obtained IC 50 values of T-1-MTA against A549 and HCT-116 malignant cells were 22.49 and 24.97 μM, respectively. T-1-MTA's anticancer potential was close to that of erlotinib.
As a confirmation of the computed safety pattern of T-1-MTA and to explore its selectivity, T-1-MTA was tested against the W138 human normal cell line. T-1-MTA showed a high IC 50 value of 55.14 μM as well as very high selectivity indexes (SI) of 2.4 and 2.2 against the two cancer cell lines, respectively (Fig 18).
2.9.3. Cell cycle analysis and apoptosis assay. Firstly, the cell cycle phases of A549 after T-1-MTA's treatment was analyzed by flow cytometry according to the reported method before [55,56]. A concentration of 22.49 μM of T-1-MTA was added to A549 cells for 72 h. Then, the cancer's cell cycle was investigated. Interestingly, T-1-MTA decreased the percentage of A549 cells in the Sub-G1 and S phases from 0.75% and 68.17% to 0.36% and 28.60%, respectively. Contraversly, in the G2/M phase, the A549 percent was significantly increased from 18.69 to 49.20 after T-1-MTA's treatment ( Table 6 and Fig 19).
To verify the apoptotic effects of T-1-MTA, the apoptosis percentage in the A549 cells was examined by Annexin V and PI double stains after it was subjected of 22.49 μM of T-1-MTA for 72 h [57,58]. Interestingly, T-1-MTA reduced the viable cancer cell count. Comparing control, T-1-MTA induced higher ratio of apoptotic cells. Also, T-1-MTA caused increased the

PLOS ONE
apoptotic cells' percentage significantly in the early stage of apoptosis (from 0.07% to 21.24%) as well as the late stage of apoptosis (from 0.73% to 37.97%). Also, the necrosis percentage was elevated to be 1.78, compared to 0.04% in the control cells (Fig 19 & Table 7). In conclusion, T-1-MTA successfully arrested the A549 cell cycle at the G2/M phase causing cytotoxic potentialities that may be connected to apoptosis.

Conclusion
According to the essential structural features of EGFR inhibitors, a new lead theobrominederived candidate, T-1-MTA has been designed. An anti-EGFR potential of the T-1-MTA was showed by molecular docking and verified by six MD simulations (over an100 ns), three MM-GBSA, and three DFT studies. Likely, computational ADMET studies indicated a general drug-likeness and safety. The biological evaluation confirmed the in silico results as T-1-MTA showed EGFR inhibitory activity with IC 50   .21g) was added to a solution of the potassium 3,7-dimethyl-3,7-dihydro-1H-purine-2,6-dione 2 (0.001 mol, 0.25g) in DMF (10 mL), and the mixture was heated in a water bath for 8 h. After being poured onto ice water (200 mL), the reaction mixture was gently stirred for certain time. To   Fig 17. In vitro EGFR-inhibition potentialities of T-1-MTA (A) and erlotinib (B).

Comp.
In vitro cytotoxicity IC 50  afford T-1-MTA (Fig 20), the obtained ppt was filtered, water washed, and crystallized from methanol.

Docking studies
Was operated for T-1-MTA by MOE2014 software. The supplementary section (S2) in S1 Data includes a detailed explanation.

MD simulations
Was operated for T-1-MTA by the CHARMM-GUI web server and GROMACS 2021 [24,59]. The supplementary section (S3) in S1 Data includes a detailed explanation.

MM-GBSA
Was operated for T-1-MTA by the Gmx_MMPBSA package [60]. The supplementary section (S4) in S1 Data includes a detailed explanation.

DFT
Was operated for T-1-MTA by Gaussian 09 and GaussSum3.0 programs. The supplementary section (S5) in S1 Data includes a detailed explanation.

ADMET studies
Was operated for T-1-MTA by Discovery Studio 4.0. The supplementary section (S6) in S1 Data includes a detailed explanation.

Toxicity studies
Was operated for T-1-MTA by Discovery Studio 4.0. The supplementary section (S7) in S1 Data includes a detailed explanation.

In vitro EGFR inhibition
Was operated for T-1-MTA by Human EGFR ELISA kit. The supplementary materials (S8) in S1 Data show a comprehensive explanation.

In vitro antiproliferative activity
Was operated for T-1-MTA by MTT procedure. The supplementary materials (S9) in S1 Data show a comprehensive explanation.

Safety assay
Was operated for T-1-MTA by MTT procedure utilizing W138 cell lines. The supplementary section (S10) in S1 Data includes a detailed explanation.

Cell cycle analysis and apoptosis
Was operated for T-1-MTA flowcytometry analysis technique. The supplementary section (S11 and S12) in S1 Data includes a detailed explanation.